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26 Abstract 

27 Global vegetation models predict rapid poleward migration of tundra and boreal forest vegetation 

28 in response to climate wanning. Local plot and air-photo studies have documented recent 

29 changes in high-latitude vegetation composition and structure, consistent with wanning trends. 

30 To bridge these two scales of inference, we analyzed a 24-year (1986-2010) Landsat time series 

31 in a latitudinal transect across the boreal forest-tundra biome boundary in northern Quebec 

32 province, Canada. This region has experienced rapid warming during both winter and summer 

33 months during the last forty years. Using a per-pixel (30 m) trend analysis, 30% of the 

34 observable (cloud-free) land area experienced a significant (p < 0.05) positive trend in the 

35 Nonnalized Difference Vegetation Index (NDVI). However, greening trends were not evenly 

36 split among cover types. Low shrub and graminoid tundra contributed preferentially to the 

37 greening trend, while forested areas were less likely to show significant trends in NDVI. These 

38 trends reflect increasing leaf area, rather than an increase in growing season length, because 

39 Landsat data were restricted to peak-summer conditions. The average NDVI trend (0.007/yr) 

40 corresponds to a leaf-area index (LAI) increase of ~0.6 based on the regional relationship 

41 between LAI and NDVI from the Moderate Resolution Spectroradiometer (MODIS). Across the 

42 entire transect, the area-averaged LAI increase was -0.2 during 1986-2010. A higher area- 

43 averaged LAI change (-0.3) within the shrub-tundra portion of the transect represents a 20-60% 

44 relative increase in LAI during the last two decades. Our Landsat-based analysis subdivides the 

45 overall high-latitude greening trend into changes in peak-summer greenness by cover type. 

46 Different responses within and among shrub, graminoid, and tree-dominated cover types in this 

47 study indicate important fine-scale heterogeneity in vegetation growth. Although our findings 

48 are consistent with community shifts in low-biomass vegetation types over multi-decadal time 
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49 scales, the response in tundra and forest ecosystems to recent warming was not unifonn. 
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50 1.0 Introduction 

51 

52 Climate exerts a primary control on the extent of forest cover and other vegetation types within 

53 Arctic and sub-Arctic ecosystems. Recent wanning has been most rapid at high latitudes, and 

54 stronger warming expected in the next century may shift the distribution of vegetation types at 

55 these latitudes (IPCC 2007). Dynamic global vegetation models (DGVMs) predict a 

56 temperature-induced growth response in high-latitude ecosystems, leading to a poleward or 

57 upslope expansion of boreal forest and an increase in the boreal forest carbon sink over the 

58 course of the 21 st century (Emanuel et al. 1985; Pastor & Post 1988; Prentice & Fung 1990; 

59 White et al. 2000; Pannesan & Yohe 2003; Lucht et al. 2006; IPCC 2007). In many of these 

60 scenarios, the vegetation response to wanning is both widespread and rapid during the 2 1 st 

61 century, which suggests that early signs of wanning-induced biome shifts might already be 

62 observable. 

63 

64 Coarse-resolution satellite data and field observations have provided intriguing evidence for 

65 climate-driven shifts in vegetation type and condition since the mid-20 th century. Vegetation 

66 index data from the coarse-resolution Advanced Very High Resolution Radiometer (AVHRR) 

67 have suggested widespread increases in high-latitude vegetation greenness and net primary 

68 productivity (NPP) since the 1980’s (Myneni et al. 1997; Goetz et al. 2005; Bunn & Goetz 2006; 

69 Pouloit et al. 2009; Wang et al. 2011). AVHRR-based studies of high-latitude greening typically 

70 use seasonally-integrated vegetation indices. Thus, it is not clear whether these satellite-based 

7 1 “greening” trends reflect increased peak-summer vegetation cover or lengthening growing 

72 seasons (including trends in spring snow cover). These studies also lack the spatial resolution to 
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73 delineate stand-level changes in vegetation composition or extent (Masek 2001) and suffer to 

74 varying degrees from issues associated with instrument calibration and sampling (e.g., Running 

75 et al. 2004; Gallo et al. 2005). 

76 

77 Field studies using plot measurements or repeat aerial photography suggest that recent climate 

78 wanning has led to expansion of shrub cover within tundra biomes (Van Wijk et al. 2004; Tape et 

79 al. 2006; Jagerbrand et al. 2006;Tremblay 2010; Ropars & Boudreau 2012). However, the extent 

80 to which such local trends contribute to or characterize larger, systemic change in Arctic and sub- 

81 Arctic vegetation remains unknown. Within the boreal forest biome, evidence for climate-driven 

82 change is less conclusive. A meta-analysis of pan-boreal treeline studies indicated that northward 

83 or altitudinal expansion of forest is evident in over half of study sites with coincident wanning 

84 trends (Harsch et al. 2009), yet variation in treeline form (e.g. diffuse, abrupt, island, and 

85 krummolz) suggests a diversity of responses to local conditions as well as climate (Harsch & 

86 Bader 2011). Recent studies demonstrate the possibility to link field observations and satellite- 

87 based trends in vegetation productivity (e.g., Beck et al. 2011); however, higher-resolution 

88 satellite observations are likely needed to directly scale individual tree or stand-scale growth 

89 responses to satellite resolution. Large area coverage with high resolution time series is also 

90 desirable since coarse or moderate resolution satellite data indicate a diversity of trends in tundra 

91 ecosystems and primarily negative (browning) trends over boreal forest in North America (Bunn 

92 & Goetz 2006; Pouliot et al. 2009; Zhao & Running 2010). 

93 

94 In this study, we use fine-scale Landsat observations to quantify vegetation changes within and 

95 among plant cover types over the past quarter-century. Our study considers a latitudinal transect 


5 



96 across the forest-tundra biome boundary in northern Quebec, a region that has experienced rapid 

97 wanning in both winter (November- April) and summer (May-October) seasons during the 

98 satellite era, from the early 1970s onwards (Fig. 1). Unlike previous satellite-based studies of 

99 vegetation greening, we carefully selected time series of peak-summer Landsat data to evaluate 

100 changes in vegetation composition and structure rather than changes in phenology. The study 

101 had three specific aims: 1) identify changes in high-latitude tundra and forest cover types over a 

102 period of pronounced warming using time series of Landsat Normalized Difference Vegetation 

103 Index (NDVI); 2) analyze the magnitude and distribution of change by cover type; and 3) assess 

104 the underlying ecological mechanisms of a trend in greenness by plant cover type and terrain 

105 attributes. 

106 
107 
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108 2.0 Materials and Methods 

109 

110 2.1 Data 

111 We assembled a Landsat time series transect across the forest-tundra biome boundary in northern 

112 Quebec to assess changes in summer vegetation cover during 1986-2010 (Fig. 1). The transect 

113 spanned 9 adjacent Landsat frames, covering an area of 260,000 km". The time series of Landsat 

114 data for each frame in the transect was selected to minimize the impact of phenology on trends in 

115 summer vegetation cover over the 24-year period (Fig. 2). We used the average growing season 

116 phenology during 2001-2006 from the MODIS phenology product (MCD12Q2, Zhang et al. 

117 2003) to select the peak growing season greenness for each Landsat data frame. Within this 

118 narrow window of peak greenness (day of year 185-215, or July 4 - August 3 in non-leap years), 

119 images were selected to minimize cloud cover and variability in the date of image acquisition. A 

120 total of 52 Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper Plus (ETM+) 

121 images were acquired from the United States Geological Survey (glovis.usgs.gov) and Canadian 

122 Centre for Remote Sensing (ccrs.nrcan.gc.ca) Landsat Data Archives, with an average of 6 

123 Landsat scenes per frame (Fig. 2). 

124 

125 Landsat data were converted from radiance to surface reflectance using the Landsat Ecosystem 

126 Disturbance Adaptive Processing System (LED APS), an automated atmospheric correction 

127 approach to account for absorption and scattering by atmospheric trace gases (O 3 , O 2 , CO 2 , NO 2 , 

128 and CH 4 ), aerosols, and water vapor (Masek et al. 2006). The LED APS system also generates a 

129 cloud mask layer. For this study, we implemented additional masking procedures for water 

130 (Band 4 reflectance <0.12), thin cirrus clouds (Band 1 reflectance >0.08), and contrails (Band 6 
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brightness temperatures). Finally, the Normalized Difference Vegetation Index (NDVI) was 
calculated from the masked red and near-infrared surface reflectance data for each scene (Tucker 
1979). Time series of NDVI for each Landsat frame were used to detect trends in forest and 
tundra vegetation during 1986-2010. 

Temperature trends in Boreal North America area were analyzed using monthly mean 
temperatures from the University of East Anglia Climatic Research Unit Time Series 3.1 (CRU 
TS3.1, http://badc.nerc.ac.uk/data/cru) . The CRU TS3.1 product is a gridded 0.5° x 0.5° product 
based on meteorological station data (see Mitchell & Jones 2005). Monthly mean temperatures 
were averaged for winter (November-April) and summer (May-October) seasons; previous 
studies have shown that warming winter temperatures may be important for recent treeline 
advances (e.g., Harsch et al., 2009). Simple linear regression was used to estimate trends in 
mean winter and summer temperatures for 1971-2008 and 1970-2009, respectively. The total 
change in winter and summer temperatures during 1970-2009 was estimated using linear trends 
on a per-pixel basis. Finally, seasonal temperature changes were averaged for each 0.5° latitude 
bin to estimate the gradients in winter and summer temperatures within the Landsat transect 
during this period. 


2.2 Trend Detection 

Trends in mid-summer NDVI were assessed on a per-pixel basis using least-squares regression. 
For a time series of n scenes, only pixels with n or n-1 observations were evaluated for trends in 
NDVI over time. Selecting pixels with one missing data value allowed for the use of some 
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154 cloud-filtered data and post-2003 Landsat ETM+ data while minimizing the potential for 

155 spurious trend detection. In regions with overlapping coverage from adjacent Landsat frames, 

156 the denser time series was selected to assess trends in NDVI over time. For each 30-meter pixel, 

157 the slope and statistical significance of the linear regression in NDVI values were evaluated 

158 using a Student’s t-test at 95% confidence level. 

159 

160 Large-scale disturbances from fire and wood harvest are common in North American boreal 

161 forests. Although climate warming may influence disturbance rates in boreal forests, the focus of 

162 this work was to detect vegetation changes in undisturbed regions. Therefore, we used two 

163 approaches to exclude large-scale disturbances from the analysis of NDVI trends. First, the 

164 earliest cloud- free Landsat MSS image for each scene (1972-1976) was used to digitize burn scar 

165 perimeters for fires prior to the start of each Landsat TM/ETM+ time series (—1986). These 

166 areas were eliminated from further analysis. Second, within the Landsat TM/ETM+ time series 

167 (1986-2010), a thresholding approach was used to eliminate areas with strong increases or 

168 decreases in NDVI (absolute changes greater than 0.08 NDVI) between successive images. 

169 

170 The sensitivity of the Landsat NDVI trend detection approach to real changes in vegetation cover 

171 depends, in part, on the uncertainty in the original radiometric observations (measurement 

172 errors). Recent improvements to Landsat calibration, including cross calibration of Landsat-4, -5, 

173 and -7 data using an absolute radiometric scale (Markham & Helder in press), have reduced 

174 uncertainties associated with comparing Landsat data from different sensors. The uncertainty of 

175 the LED APS atmospherically corrected products is the greater of 0.5% absolute reflectance or 

176 5% of the recorded reflectance value (la), similar to the NASA Moderate Resolution Imaging 
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177 Spectroradiometer (MODIS) sensors (Masek et al. 2006). The resulting uncertainty in any given 

178 NDVI observation is thus -0.02 (la). Based on Monte Carlo modeling of MODIS NDVI trend 

179 detection using noisy time series (Wang et al. 2012), we expect that actual NDVI trends greater 

180 than ±0.003 NDVI yr' 1 can be reliably mapped given the Landsat measurement errors and the 

181 typical number of observations for each frame in the transect. 

182 

183 Increases in NDVI over time are generally associated with increases in leaf area index (LAI) 

184 (Turner et al. 1999). Because changes in LAI may be more easily compared to modeled or 

185 measured vegetation changes, we developed relationships between NDVI and LAI using the 

186 MODIS/Aqua LAI product (MYD15A2, Knyazikhin et al. 1999) and NDVI product 

187 (MYD13A1, Huete et al. 2002). First, trends in 2002-2010 MODIS LAI and NDVI were 

188 assessed independently across the northern United States and Canada, using the same linear 

189 regression and t-test approach as the analysis of trends in Landsat NDVI. Data quality layers 

190 were used to restrict the analysis to LAI estimates from the radiative transfer algorithm and 

191 highest quality NDVI observations. Pixels with statistically significant linear trends in both 

192 MODIS LAI and MODIS NDVI were selected to derive estimates of the change in LAI per unit 

193 NDVI. The larger (continental) geographic area allowed over 200,000 MODIS NDVI-LAI pairs 

194 to be collected for comparison. 

195 

196 2.3 Spatial Analysis of NDVI Trends 

197 Trends in Landsat NDVI during 1986-2010 were assessed by cover type using forest and tundra 

198 land cover classifications derived from circa 2000 Landsat imagery. North of the tree line, 

199 designated as the Riviere aux Feuilles, tundra cover types were classified based on the CCRS 
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200 Northern Land Cover of Canada dataset 

201 (http://www.ccrs.nrcan.gc.ca/optical/landcover2000_e.php). The Canadian Forest Service's 

202 Earth Observation for Sustainable Development of Forests (EOSD) classification was used to 

203 analyze trends in NDVI by cover type for portions of the time series transect south of tree line 

204 (http://www4.saforah.org/eosdlcp/nts_prov.html). Land cover information from the two 

205 classification products was merged to create a harmonized classification for the study region with 

206 six vegetation classes, barren or exposed bedrock, and water (Table 1). Trends were further 

207 analyzed by latitude by comparing the mean NDVI trend between adjacent (north-south) scenes. 

208 Trends in NDVI and temperature were tested for correlation with latitude using Pearson’s 

209 product-moment correlation coefficient. 

210 

211 Finally, trends in Landsat NDVI for each cover type were further evaluated by slope, aspect, and 

212 elevation. Topographic infonnation was derived from the Shuttle Radar Topography Mission 

213 (SRTM) 3 Arc Second Filled Finished-B product (USGS 2006; www.landcover.org). The 

214 relationship between slope, elevation, and aspect with positive non-disturbed NDVI trend 

215 occurrence was explored for a random sample of positive and no-trend observations (equal to 

216 10% of total observations of each) using a binomial generalized linear model (R version 2.11.0, 

217 R Core Development Team 2010). 

218 

219 

220 
221 
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222 3. Results 

223 The Landsat study area covered 26 million ha in a transect across the forest-tundra biome 

224 boundary. Approximately 70% of the transect had a sufficient number of cloud-free Landsat 

225 observations to assess trends in peak-season NDVI. Of this “observable” portion, one-third 

226 (34%) of the area experienced a statistically significant trend in NDVI during 1986 - 2010 based 

227 on the T-test criterion (Fig. 3a). The remaining observable area either had small NDVI trends 

228 that were not significantly different from zero or exhibited significant year-to-year variability 

229 that precluded statistical confidence. 

230 

231 Almost all of the statistically significant NDVI trends were positive. Large-scale forest 

232 disturbance events prior to or during the study period, including forest fires and timber harvests, 

233 accounted for 3.2% of the area with significant NDVI trends. Excluding these disturbances, 

234 98.95% of the remaining trend area had positive (greening) trends and only 1.05% of the trend 

235 area exhibited negative (browning) NDVI trends (Fig. 3a). The mean positive trend was an 

236 increase 0.007 NDVI yr' 1 for a total increase of 0. 17 NDVI over the entire 24-year time series 

237 (Fig. 3b). Positive NDVI trends were concentrated north of the treeline; nearly half (48%) of the 

238 observed area north of treeline had a statistically significant positive NDVI trend compared to 

239 only 25% of observable area south of treeline. Latitude and the frequency of statistically 

240 significant NDVI trends were positively correlated (Pearson’s R = 0.819, P = 0.0248), although 

241 the magnitude of the NDVI trend was not correlated with latitude (Fig. 4, Table 2). Temperature 

242 trends, both annual and cumulative, were not significantly correlated with trends in NDVI 

243 abundance or magnitude (Table 2). 

244 
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245 NDVI trends showed significant associations with specific land cover types (Fig. 5, Table 3). 

246 Low shrubs and graminoid tundra contributed preferentially to the observed greening trend. Out 

247 of the area showing a positive NDVI trend, 38% occurred in regions classified as either 

248 low/dwarf shrubs or graminoid tundra, even though these types only make up 22% of the study 

249 area. Within the total area of tall shrubs, 30% of observations showed a positive trend. The rate 

250 of NDVI increase for all cover types was broadly similar, varying from 0.0055 NDVI year _1 for 

251 sparsely vegetated areas to 0.0075 NDVI yr' 1 for graminoid tundra. Although forests comprise 

252 the majority of the southern portion of the transect, they contributed less than 10% to the overall 

253 greening trend. Within forests, 15% of total observations showed a positive trend, with a mean 

254 increase of 0.0064 NDVI yr' 1 . In contrast, 50-60% of the low/dwarf shrub and graminoid tundra 

255 areas showed significant NDVI increases. 

256 

257 Positive trends in peak summer NDVI correspond to an increase in LAI over time. Using 

258 MODIS data for the study region, statistically significant changes in NDVI between 0.005 - 0.01 

259 yr' 1 corresponded to changes in LAI of -0.02 (mode) to 0.03 (median) LAI yr' 1 (Fig. 6). This 

260 relationship between MODIS LAI and NDVI was invariant across the multiple cover types 

261 occurring in the study region. We applied a value of 0.025 LAI yr' 1 to those Landsat NDVI 

262 trends not associated with disturbance, assuming that areas not exhibiting significant trends 

263 experienced no change in LAI. The resulting estimate of the area-averaged LAI change across 

264 the transect was 0.2 LAI over the 24-year record. However, for tundra cover types in the 

265 northern portion of the transect, the estimated 24-year LAI increase was ca. 0.3. Tundra and 

266 shrub-tundra LAI values generally fall in the range of 0.5 to 1.5 (Asner et al. 2003; Beringer et 
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267 al. 2005; Williams 2008), thus an LAI increase of 0.3 translates to roughly a -20-60% increase 

268 in leaf area over two decades. 

269 Geographically, the highest magnitude of NDVI change occurred proximal to and north of the 

270 regional treeline, which roughly coincides with the Riviere aux Feuilles (Fig. 3b). Greening 

271 trends were less abundant in the forested regions in the southern half of the transect, and trends 

272 decayed in frequency and magnitude along the northern edge bordering the Hudson Strait. 

273 

274 Topography was an important predictor of greening trends over the study domain (Table 4). 

275 While the landscape is generally flat, locations with higher slopes and elevations were negatively 

276 correlated with the frequency of detecting greening trends. North- and northeast- facing slopes 

277 were least likely to exhibit a positive trend, and western and southwestern facing slopes were the 

278 most likely. Stronger NDVI trends were detected along two major rivers, the Riviere aux 

279 Feuilles and the southern reaches of the Riviere Amaud. Topographic associations between 

280 valley bottoms and vegetation growth likely reflect more favorable edaphic conditions along the 

281 channel banks, as well as more sheltered microclimates and available water. 

282 

283 Most cover types exhibited consistent linear changes in NDVI over the 24-year study period. 

284 The temporal distribution of the Landsat images did not support a regional, year-by-year analysis 

285 of greening trends. Instead, we divided observations into “early” (1985-1990), “mid” (1998- 

286 2001), and “late” (2008-2010) intervals to evaluate rates of NDVI change by cover type over 

287 time (Fig. 7). Most classes maintained similar rates of NDVI change during both early and late 

288 intervals. However, in wetland and tall shrub classes, NDVI increases were slower in the later 

289 period. 

290 
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291 4. Discussion 

292 

293 Using time series of Landsat data, we found a strong mid-summer greening trend across the 

294 northern Quebec area. This trend corresponds to significant increases in peak growing season 

295 leaf area. Graminoid and shrub-tundra classes contributed nearly 60% of the greening trends 

296 identified in this study. These low-biomass vegetation types experienced a 20-60% relative 

297 increase in green leaf area over the 24-year study period — a rapid and significant increase 

298 relative to existing phytomass. Whether these LAI gains reflect the growth of individual plants, 

299 an increase in the density of individuals, or community-scale changes such as shifts from 

300 graminoid- to shrub-dominated tundra remains an important area for further study. Large and 

301 persistent changes in LAI identified in this study may also alter biophysical feedbacks, including 

302 seasonal changes in albedo (e.g., Randerson et al. 2006; Bonan 2008), snow cover, and turbulent 

303 fluxes (e.g., Lee et al. 2011). 

304 

305 Changes in shrub and graminoid tundra in this study are consistent with regional evidence of 

306 shrub expansion in tundra ecosystems (Ropars & Boudreau, 2012) and similar findings across 

307 North America (Chapin 1995, Sturm et al. 2001; Tape et al. 2006, Jagerbrand 2005; Van Wijk et 

308 al. 2004). A recent air photo analysis for a region west of our transect revealed increases in 

309 dwarf birch ( Betula glandulosa Michx.) shrub cover up to 47% over a 5 1-year period (1957- 

310 2008), with larger changes on low-altitude sandy terraces than on exposed hilltops (Ropars & 

311 Boudreau, 2012). The Landsat data suggest that this shrub expansion is not isolated, but 

312 widespread across the tundra region of northern Quebec. Quebec and Alaska were the only 

313 regions in North America with strong wanning trends in both winter and summer semesters 
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during 1970-2009, and both regions have consistent reports of increasing shrub cover during the 
satellite era. We did not identify strong correlations between the magnitude of recent 
temperature changes and increases in fractional cover, possibly due to consistent and strong 
(1.5°-2.5°) warming across the entire study region. Given favorable climatic conditions, 
landscape heterogeneity and species-level responses may be stronger predictors of vegetation 
change. 

Biophysical mechanisms operating at the local scale may contribute to observed ND VI/LAI 
increases in tundra cover types. First, snowdrifts are more likely to form at lower elevations, 
trapping organic debris and leaf litter (Fahnestock et al. 2000). Snow protects newly-established 
shrubs from harsh winter conditions, and wanner soil temperatures under deep snow may 
increase microbial activity that mobilizes additional nutrients for shrub growth (Sturm et al. 
2005). A positive feedback between shrub cover and snow is consistent with larger LAI 
increases at lower slopes and elevations in this study. In more sparsely vegetated cover types, 
dominated by lichen, mosses, and bryoids, another feedback cycle between the lichen-caribou- 
woody plant communities may be important. Caribou trampling destroys the lichen and exposes 
mineral soil, allowing for an increase in seedling establishment for dominant subarctic tree and 
shrub species such as black spruce ( Picea mariana; Dufour-Tremblay & Boudreau 2011) and 
dwarf birch ( B . glandulosa; Ropars & Boudreau, 2012). The activity of the Leaf River Caribou 
Herd in western subarctic Quebec peaked during the mid-1990s to mid-2000s (Dufour-Tremblay 
& Boudreau, 2011; Alexandre Truchon-Savard, pers. comm.), suggesting that browsing and soil 
disturbances from large herbivores may contribute to the patterns of shrub and graminoid tundra 
changes in this study. The availability of bare soil, whether from caribou disturbance or other 


16 


337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 


disturbances like frost boils or cryoturbation, combined with milder conditions more favorable 
for seed production and seedling establishment, may allow for the encroachment of woody 
species and other vascular vegetation into sparsely vegetated areas. 

This study found less conclusive evidence for vegetation changes within forest areas. While 
most observations of NDVI trend within the forested parts of the transect were positive, a much 
smaller area showed a statistically significant trend compared to graminoid- and shrub- 
dominated regions. In contrast to the shrub expansion studies, evidence for northern advance of 
treeline into tundra has been mixed (Harsch et al. 2009). At the treeline, a positive temperature 
trend may not necessarily correlate with the northward expansion of trees, given the influence of 
water availability, soil properties, competition, or pests on the spatial arrangement of trees 
(Meunier et al. 2007). Furthermore, the responsiveness of non-tree species in forest 
communities, such as shrubs, to a positive temperature trend may be suppressed by tree cover 
(Boudreau & Villenueve-Simard 2012). Although within-stand changes in forest leaf area were 
less common, it is possible that expansion of tree species into tundra communities dominated by 
tall shrubs or other functional groups contributed some of the observed changes in other 
vegetation classes reported in the study. Additional field studies in areas of recent change are 
needed to identify the contributions from the growth of existing individuals (Gamache & Payette 
2004; Beck et al. 2011) and establishment of new individuals (e.g., Danby & Hik 2007) to recent 
increases in vegetation cover. 

Observed NDVI trends in forest may reflect recovery from historic disturbances, despite efforts 
to mask out large-scale disturbances from fires and forestry operations visible during the Landsat 
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360 era. In eastern Canada, severe fires sharply reduce LAI, and vegetation regrowth occurs over 

361 century timescales, during which post- fire succession is likely to overshadow climate-driven 

362 trends in vegetation (Girard et al. 2008). Two other factors may contribute to the lack of 

363 observed forest greening within a time series of Landsat data. First, the establishment and 

364 growth of trees is an inherently slower process compared to growth of existing individuals (e.g., 

365 Danby & Hik 2007), and the 24-year Landsat record may not be long enough to identify changes 

366 within forest stands. Second, forested areas tend to have higher initial NDVI values. Since 

367 NDVI saturates at modest LAI values (~3.0), small LAI increases within existing forest stands 

368 may not be obvious from the remote sensing data. 

369 

370 DGVMs suggest poleward migration of biomes as a long-term response to climate warming 

371 (Lucht et al. 2006). The observed association between shrub cover types and increased NDVI is 

372 generally consistent with the concept that woody plants can take advantage of warmer conditions 

373 and grow more vigorously. In areas of mixed graminoid and shrub cover types, the competitive 

374 advantage of shrubs should lead to a long-term shift in composition, and ultimately a poleward 

375 shift in the biome boundary. However, the satellite data do not yet provide unambiguous 

376 evidence for geographic biome shifts as opposed to simply increasing LAI within existing biome 

377 distributions. As noted by others, the ability of vegetation communities to expand their range 

378 depends not just on increased productivity, but on overcoming a host of ecological constraints 

379 (Rozensweig et al. 2008). Particularly in the boreal environment of Canada, small lakes and 

380 rocky outcrops present innumerable fine-scale barriers to propagation and expansion. The poor 

381 reproductive capacity of frontier tree species such as Picea mariana may also somewhat explain 

382 observed lags between warming and vegetation growth, both above and below the subarctic 
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treeline (Gamache & Payette 2004), although an increase in seed viability was noticed near the 
treeline in recent years (Dufour-Tremblay & Boudreau, 2011). These fine-scale barriers to forest 
expansion constitute macro-scale “resistance” to biome shifts that are not considered in the 
current generation of DGVMs. 

Our results complement previous studies of high-latitude vegetation change using moderate or 
coarse-resolution satellite data (e.g., Pouliot et al. 2009). We used Landsat time series to 
subdivide the overall greening trend into increases in LAI for specific cover types. The Landsat- 
based approach in this study could be expanded to evaluate climate-driven shifts in vegetation in 
other regions, within the limits of the existing Landsat data archive (Goward et al. 2006). High- 
resolution time series over the 35+ year Landsat record provide invaluable observational data to 
refine and benchmark ecological models. However, there are several important limitations of 
this work that could be addressed in future studies. First, given the uncertainties in the MODIS 
LAI product in high-latitude areas, the uncertainties in the derived MODIS NDVI-LAI 
relationship, and the difficulties of scaling field-observed LAI to satellite resolution, area- 
averaged LAI increases in this study should be interpreted cautiously. Second, we evaluated 
trends in Landsat NDVI by cover type using vegetation classification data from a snapshot in 
time (circa 2000). Classification information from 2000 may already incorporate growth of 
woody vegetation during 1986-2000, such that areas classified as shrublands in 2000 had lower 
amounts of woody cover at the beginning of the study period. Third, multispectral remote 
sensing has limited sensitivity to subtle changes in composition and structure, especially in 
closed-canopy forest conditions. The addition of hyperspectral imagery (to map compositional 
gradients) and LiDAR(to map structure) would provide a more comprehensive benchmark of 
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406 current conditions for future studies of climate-driven vegetation changes. Finally, we identified 

407 linear changes in NDVI over the Landsat study domain. Non-linear greening or browning 

408 responses to recent climate warming may therefore be underestimated in this study. 

409 

410 Using time series of Landsat observations, we mapped widespread vegetation greening in 

411 northern Quebec over the last 24 years. The observed NDVI increases were concentrated in 

412 graminoid and shrub-tundra areas, leading to an area-averaged LAI increase of ~0.2 across the 

413 entire transect, or -0.3 for the northern tundra-dominated portion. The latter figure represents a 

414 20-60% relative increase compared to typical shrub-tundra LAI values. These findings expand 

415 the spatial extent of previous field and air photo studies used to characterize changes in shrub 

416 cover. Our results also provide a fine-scale evaluation of the contribution of different cover 

417 types to trends detected from coarse-resolution satellite data. The coincidence of the shrub 

418 greening trend with an area of rapid winter and summer wanning supports the hypothesis that 

419 wanner temperatures favor the growth of woody plants at high latitudes (Sturm et al. 200 1 ; 

420 2005; Tape et al. 2006). In contrast, positive NDVI trends within forested areas were less 

421 common, suggesting that the forest response to recent wanning may be occuning more slowly, 

422 or that Landsat data alone may be insufficient to identify growth responses in these ecosystems 

423 and additional data (e.g., LiDAR) may be needed to characterize temperature-induced vegetation 

424 changes within boreal forest communities. 

425 

426 


20 



427 

428 

429 

430 

431 

432 

433 

434 

435 

436 


6. Acknowledgements 

This work was supported by the NASA Terrestrial Ecology Program. We thank M.A. Wulder for 
guidance on land cover information for eastern Canada, and we appreciate comments and 
suggestions on a previous version of this manuscript from H. Margolis and an anonymous 
reviewer. 


21 



437 7. References 

438 

439 Asner GP, Schurlock J, and Hicke JA (2003) Global synthesis of leaf area index observations: 

440 implications for ecological and remote sensing studies. Global Ecology & Biogeography, 12, 

441 191-205. 

442 

443 Beck SA, Homing S, Goetz SJ, Loranty MM, Tape KD (2001) Shrub Cover on the North Slope 

444 of Alaska: a circa 2000 Baseline Map. Arctic, Antarctic, and Alpine Research 43, 355-363. 

445 

446 Beringer J, Chapin FS, Thompson CC, McGuire AD (2005) Surface energy exchanges along a 

447 tundra-forest transition and feedbacks to climate. Agricultural and Forest Meteorology, 131, 

448 143-161. 

449 

450 Bonan GB (2008) Forests and Climate Change: Forcings, Feedbacks, and the Climate Benefits of 

451 Forests. Science, 320, 1444-1449. 

452 

453 Boudreau S & Villeneuve-Simard MP (2012) Dendrochronological evidence of shrub growth 

454 suppression by trees in a subarctic lichen woodland. Botany 90, 151-156. 

455 

456 Bunn AG & Goetz SJ (2006) Trends in satellite-observed circumpolar photosynthetic activity 

457 from 1982 to 2003: the influence of seasonality, cover type, and vegetation density. Earth 

458 Interactions, 10, 1-19. 

459 

460 Chapin FS III, Shaver GR, Giblin AE, Nadelhoffer KJ, and Laundre JA (1995) Responses of 

461 arctic tundra to experimental and observed changes in climate. Ecology, 76, 694-7 1 1 . 

462 

463 Danby RK & Hik DS (2007) Variability, contingency and rapid change in recent subarctic alpine 

464 tree line dynamics. Journal of Ecology, 95, 352-363. 

465 

466 Doufour-Tremblay G & Boudreau SB (2011) Black spruce regeneration at the treeline ecotone: 

467 synergistic impacts of climate change and caribou activity. Canadian Journal of Forest, 41, 460- 

468 468. 

469 

470 Emanuel WR, Shugart HH, Stevenson MP (1985) Climate change and the broad scale 

471 distribution of terrestrial ecosystem complexes. Climatic Chang, e 7: 29-43. 

472 

473 Gallo K, Ji L, Reed B, Eidenshink J, Dwyer J (2005). Remote Sensing of the Environment, 99, 

474 221-231. 

475 

476 Gamache I & Payette S (2004) Height growth response of tree line black spruce to recent climate 

477 warming across the forest-tundra of eastern Canada. Journal of Ecology. 92, 835-845. 

478 

479 Girard F, Payette S, Gagnon R (2008) Rapid expansion of lichen woodlands within the closed- 

480 crown boreal forest zone over the last 50 years caused by stand disturbances in eastern Canada. 

481 Journal of Biogeography. 35, 529-537. 

482 


22 



483 

484 

485 

486 

487 

488 

489 

490 

491 

492 

493 

494 

495 

496 

497 

498 

499 

500 

501 

502 

503 

504 

505 

506 

507 

508 

509 

510 

511 

512 

513 

514 

515 

516 

517 

518 

519 

520 

521 

522 

523 

524 

525 

526 

527 

528 


Goetz SJ, Bunn AG, Fiske GJ, Houghton RA (2005) Satellite-observed photosynthetic trends 
across boreal North America associated with climate and fire disturbance. Proceedings of the 
National Academies of Science, 102 , 13521-5. 

Goward S, Arvidson T, Williams D, Faundeen J, Irons J, Franks S (2006) Historical record of 
Landsat global coverage: Mission operations, NSLRSDA, and international cooperator stations. 
Photogrammetric Engineering & Remote Sensing, 72 , 1155-1169. 

Harsch MA, Hulme PE, McGlone MS, Duncan RP (2009) Are treelines advancing? A global 
meta-analysis of treeline response to climate warming. Ecology Letters, 12 , 1040-1049. 

Harsch MA & Bader MY (2011) Treeline form - a potential key to understanding treeline 
dynamics. Global Ecology and Biogeography, 20:582-596. 

Huete A, Didan K, Miura T, Rodriguez EP, Gao X, Ferreira LG (2002) Overview of the 
radiometric and biophysical performance of the MODIS vegetation indices. 195-213. 

IPCC (2007) Climate change 2007:Impacts, adaptation and vulnerability. Contribution of 
Working Group II to the Fourth Assessment Report of the Intergovernmental Panel on Climate 
Change Parry ML, Canziani OF, Palutikof JP, van der Linden PJ, Hanson CE, eds. Cambridge 
University Press: 211-272. 

Jagerbrand AK, Lindblad KEM, Bjork RG, Alatao JM, Molau U (2006) Bryophyte and lichen 
diversity under simulated environmental change compared with observed variation in 
unmanipulated alpine tundra. Biodiversity and Conservation, 15 , 4453-4475. 

Knyazikhin Y, Glassy J, Privette JL, et al. MODIS Leaf Area Index (LAI) and Fraction of 
Photosynthetically Active Radiation Absorbed by Vegetation (FPAR) Product (MOD 15) 
Algorithm, Theoretical Basis Document, 126 pp., NASA Goddard Space Flight Center, 
Greenbelt, Md., 1999. 

Kullman L (2001) 20th Century Climate Warming and Tree-Limit Rise in the Southern Scandes 
of Sweden. Ambio, 30 , 72-80. 

Lee X, Goulden ML, Hollinger DY, et al. (2011) Observed increase in local cooling effect of 
deforestation at higher latitudes. Nature, 479 , 384-387. 

Lucht W, Schaphoff S, Erbrecht T, Heyder U, Cramer W (2006). Terrestrial vegetation 
redistribution and carbon balance under climate change. Carbon Balance and Management, 1 , 
1 - 6 . 


Markham BL & Helder D (in press) Forty-year calibrated record of earth-reflected radiance from 
Landsat: a review, Remote Sensing of the Environment. 

Masek JG (2001) Stability of boreal forest stands during recent climate change: evidence from 
Landsat satellite imagery. Journal of Biogeography, 28 , 967-976. 


23 


529 

530 Masek JG, Vermote EF, Saleous N, et al. (2006) A Landsat surface reflectance data set for North 

531 America, 1990-2000. Geoscience and Remote Sensing Letters, 3, 68-72. 

532 

533 Meunier C, Sirois L, Begin Y (2007) Climate and Picea mariana seed maturation relationships: a 

534 multi-scale perspective. Ecological Monographs, 77, 361-376. 

535 

536 Mitchell TD & Jones PD (2005) An improved method of constructing a database of monthly 

537 climate observations and associated high-resolution grids. International Journal of Climatology, 

538 25,693-712. 

539 

540 Myneni RB, Keeling CD, Tucker CJ, Asrar G (1997) Increased plant growth in the northern high 

541 latitudes from 198 lto 1991. Nature, 386, 698-702. 

542 

543 Pastor J & Post WM (1988) Response of northern forests to CCh-induced climate change. 

544 Nature, 334 , 55-58. 

545 

546 Parmesan C & Yohe G (2003) A globally coherent fingerprint of climate change impacts across 

547 natural systems. Nature, 421 , 37-42. 

548 

549 Pouliot D, Latifovic R, Olthof I (2009) Trends in vegetation NDVI from 1 km AVHRR data over 

550 Canada for the period 1985-2006. International Journal of Remote Sensing, 30 , 149-168. 

551 

552 Prentice KC & Fung IY (1990) The sensitivity of terrestrial carbon storage to climate change. 

553 Nature, 346 , 48-50. 

554 

555 Randerson JT, Fiu H, Flanner MG, et al. (2006) The impact of boreal forest fire on climate 

556 warming. Science, 314 , 1130-1132. 

557 

558 Ropars P, & Boudreau S (2012) Shrub expansion at the forest-tundra ecotone: spatial 

559 heterogeneity linked to local topography. Environmental Research Letters 7 (015501). 

560 

561 Rosenzweig C, Karoly D, Vicarelli M, et al. (2008) Attributing physical and biological impacts 

562 to anthropogenic climate change. Nature, 453 , 353-357. 

563 

564 Running SW, Nemani RR, Heinsch FA, Zhao M, Reeves M, Hashimoto H (2004). A continuous 

565 satellite-derived measure of global terrestrial primary production. BioScience, 54 , 547:560. 

566 

567 Sturm M, Racine C, Tape K (2001) Climate change - increasing shrub abundance in the Arctic. 

568 Nature, 411 , 546-547. 

569 

570 Sturm M, Schimel J, Michaelson G, Welke JM (2005) Winter biological processes could help 

571 convert Arctic tundra to shrubland. Bioscience, 55,17-26. 

572 

573 Tape K, Sturm M, Racine C (2006) The evidence for shrub expansion in Northern Alaska and the 

574 Pan-Arctic. Global Change Biology, 12 , 686-702. 


24 



575 

576 

577 

578 

579 

580 

581 

582 

583 

584 

585 

586 

587 

588 

589 

590 

591 

592 

593 

594 

595 

596 

597 

598 

599 

600 

601 

602 

603 

604 

605 

606 

607 

608 

609 

610 

611 

612 

613 

614 

615 


Tremblay B (2010) Augmentation recente du couvert ligneux erige dans les environs de 
Kangiqsualujjuaq (Nunavik, Quebec). MSc thesis, Universite du Quebec a Trois-Rivieres, Trois- 
Rivieres, Quebec, Canada. 

Tucker CJ (1979) Red and Photographic Infrared Linear Combinations for Monitoring 
Vegetation, Remote Sensing of Environment, 8,127-150. 

Turner DP, Cohen WB, Kennedy RE, Fassnacht KS, & Briggs JM (1999) Relationships between 
leaf area index and Landsat TM spectral vegetation indices across three temperate zone sites. 
Remote Sensing of Environment, 70 , 52-68. 

University of East Anglia Climatic Research Unit (CRU) [Phil Jones, Ian Harris]. CRU Time 
Series (TS) high resolution gridded datasets, [Internet]. NCAS British Atmospheric Data Centre, 
2008. http://badc.nerc.ac.uk/view/badc.nerc.ac.uk ATOM dataent 1256223773328276 

Van Wijk MT, Clemmensen KE, Shaver GR, et al. (2004) Long-term ecosystem level 
experiments at Toolik Lake, Alaska, and at Abisko, Northern Sweden: generalizations and 
differences in ecosystem and plant type responses to global change. Global Change Biology, 10 , 
105 -123 . 

Wang D, Morton D, Masek J, et al. (2012) Impact of sensor degradation on the MODIS NDVI 
time series. Remote Sensing of Environment. 119 , 55-61. 

Wang X, Piao S, Ciais P, Li J, Friedlingstein P, Koven C, Chen A (2011) Spring temperature 
change and its implication in the change of vegetation growth in North America from 1982 to 
2006. Proceedings of the National Academies of Science, 108 , 1240-1245. 

White A, Cannell MGR, Friend AD (2000) The high-latitude terrestrial carbon sink: a model 
analysis. Global Change Biology, 6 , 227-245. 

Williams M, Bell R, Spadavecchia L, Street LE, van Wijk MT (2008) Upscaling leaf area index 
in an Arctic landscape through multi-scale observations. Global Change Biology, 14 , 1517-1530. 

Zhang X, Friedl MA, Schaaf CB, Strahler AH, Hodges JCF, Gao F, Reed B, Huete A (2003) 
Monitoring vegetation phenology using MODIS. Remote Sensing of Environment , 84 , 471-475. 

Zhao M and Running SW (2010) Drought -Induced Reduction in Global Terrestrial Net Primary 
Production from 2000 Through 2009. Science, 329 , 940 - 943. 


25 


616 


617 

618 

619 

620 
621 
622 

623 

624 

625 

626 
627 


Table 1. Harmonization scheme for a simplified vegetation classification of the study area, 
based upon the Canadian Centre for Remote Sensing (CCRS) Northern Land Cover of 
Canada and the Canadian Forest Service Earth Observation for Sustainable Development of 
Forests (ESOD) datasets. 


CCRS classification 

EOSD classification 

Harmonized Classification 


Herbaceous (grasses, crops, forbs, graminoids; 

Graminoid 

Tussock graminoid tundra 
Wet sedge 

20% ground cover) 


Moist to dry non-tussock graminoid/ dwarf 
shrub tundra 

Prostrate dwarf shrub 

Low shrub (< 40cm; >25% cover) 


Low & dwarf shrub 

Tall shrub (>40cm; >25% cover) 

Shrub- tall (>2m; 20% ground cover) 
Shrub- low (<2m; 20% ground cover) 

Tall shrub 


Coniferous: Dense, Open, and Sparse 1 
Broadleaf: Dense, Open, and Sparse 1 
Mixed Wood: Dense, Open, and Sparse 1 

Forest 


Bryoids (bryophytes and lichen; 20% ground 

Sparse vegetation 

Sparsely vegetated bedrock 

Sparsely vegetated till-colluvium 

Bare soil with cryptogam crust-frost boils 

cover or 1/3 of vegetation) 


Wetlands 

Wetland- Coniferous 
Wetland- Broadleaf 
Wetland- Mixed Wood 
Wetland- Shrub- Tall 
Wetland- Shrub- Low 

Wetlands 

Barren 

Rock/rubble 

Exposed land (<5% vegetation) 

Barren & exposed surfaces 

Water 

Water 

Water 

Ice/snow 

Cloud 

No data 

Shadow 

Shadow 


No data 

Snow/ice 
No data 



1 Coniferous, broadleaf, and mixed wood classes were further subdivided in the EOSD classification as dense (>60% 
crown closure), open (26-60% crown closure), and sparse (10-25% crown closure). 
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Table 2. Correlation coefficients for the relationships among latitude, summer 
and winter temperature (T) changes (1970-2009), and Landsat NDV1 trends 
(1986-2010), summarized at 0.5° resolution. 


Variables 

Latitude 

Winter T Change 

Summer T Change 

Latitude 

- 



Winter T Change 

0.39 

- 

- 

Summer T Change 

0.47 

- 

- 

Fraction + NDVI Trend 

0.82* 

0.19 

0.28 

Fraction - NDVI Trend 

-0.71 

-0.56 

-0.64 

Fraction no NDVI Trend 

-0.81* 

-0.17 

-0.26 

Positive NDVI Trend 

-0.62 

-0.04 

-0.19 

Negative NDVI Trend 

-0.34 

-0.23 

-0.38 

* p <0.05; **p <0.01; ***p 

<0.001 
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Table 3. Positive Landsat NDV1 trends by land cover class for the Quebec study region. 


Class 

Fraction of 
study area 


Fraction of positive 
NDVI trend area 

Fraction of class with 
positive NDVI trend 

Mean annual positive 
NDVI trend (±1S.D.) 

Barren and exposed 
surfaces 


8.0 

8.1 

35.0 

6.0 x 10" 2 ± 2.2 x 10" 2 

Sparse vegetation 


16.5 

15.7 

32.7 

5.5 x 10" 2 ± 2.8 x 10" 2 

Tall shrubs 


25.6 

23.1 

31.0 

7.5 x 10" 2 ± 3.5 x 10" 2 

Wetlands 


2.5 

2.1 

28.7 

5.9 x 10" 2 ± 3.9 x 10" 2 

Forest 


20.9 

9.7 

16.0 

6.4 x 10" 2 ± 3.8 x 10' 2 

Low and dwarf shrubs 


14.1 

23.1 

56.5 

6.3 x 10" 2 ± 2.3 x 10' 2 

Graminoid 


8.3 

14.6 

60.8 

7.5 x 10" 2 ± 2.7 x 10’ 2 

Water 


4.2 

3.5 

28.2 

5.6 x 10" 2 ± 2.4 x 10’ 2 
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Table 4. Topographic coefficients derived from 
generalized linear model (GLM) fit of greening trends. 
Based upon Landsat data from 1985-2010 for northern 
Quebec. 

Binomial GLM: Positive 
non-disturbed trend/no trend. 



Estimate c 

Intercept 

1.282 

Slope 3 

-0.03651 

Elevation 3 

-0.006545 

Aspect b : 


Flat 

-0.5449 

N 

-0.08753 

NE 

-0.1102 

E 

-0.1271 

SE 

-0.07817 

S 

-0.04362 

sw 

0.05415 

w 

0.08525 

NW 

NA 

a continuous; b ordinal 


c p-values for all variables < 2e-16 
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649 

650 Fig. 1. Changes in mean winter (left) and summer (right) temperatures between 1970 and 2009 

651 across Boreal North America based on the CRU TS3.1 dataset. The locations of the Landsat 

652 transect (white) and boreal forest biome (green) are also shown. 

653 
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654 Figure 2. The temporal distribution of the Landsat data (1985-2010) used in this study. Each 

655 panel represents a frame (path-row location), and the images used are shown by year (x-axis) and 

656 day of year (y-axis). Images were selected within a window of peak greenness (day of year 185- 

657 215) whenever high quality, minimal cloud-covered images were available. Closed circles 

658 indicate Landsat-5 TM data; open circles indicate Landsat-7 ETM+ data. 

659 

660 Figure 3. a. Locations of positive, negative, and no NDVI trend across the study region, based 

661 upon Landsat data from 1985-2010 for northern Quebec. White regions signify that data were not 

662 in sufficient quantity to detennine a statistical trend. Red regions denote areas of known 

663 disturbances; b. The magnitude of trend across the study region. 

664 

665 Figure 4. The mean annual NDVI trend by latitude (top). Trends in winter and summer 

666 temperatures (1970-2009) and fraction of NDVI change by latitude (bottom). 

667 

668 

669 Figure 5. (left) The fraction of observable (non cloud) study area experiencing positive, 

670 negative, or no trend in Landsat NDVI (1986-2010); (right) distribution of land cover types 

67 1 among the area experiencing positive (greening) trend. 

672 

673 Figure 6. (a) Relationship between MODIS and Landsat NDVI (aggregated to 500m resolution) 

674 for pixels showing a statistically significant, positive NDVI trends; (b) Relationship between 

675 MODIS NDVI and LAI trends for the northern US and Canada. For the range of NDVI changes 

676 considered in this study (0.005-0.01 NDVI/yr), the corresponding modal and median values of 

677 LAI change are 0.02 and 0.03 LAI/yr, respectively. 

678 

679 Figure 7. Mean NDVI values per class for 1986, 2000, and 2010 reference periods. Data for 

680 each Landsat frame in the time series transect were selected within ±2 years of these reference 

681 years. 
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